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We present a wavefunction-based approach to correlated ab 
initio calculations on crystalline insulators of infinite extent. 
It uses the representation of the occupied and the unoccupied 
(virtual) single-particle states of the infinite solid in terms of 
Wannier functions. Electron correlation effects are evaluated 
by considering virtual excitations from a small region in and 
around the reference cell, keeping the electrons of the rest 
of the infinite crystal frozen at the Hartree-Fock level. The 
method is applied to study the ground state properties of the 
LiH crystal, and is shown to yield rapidly convergent results. 



I. INTRODUCTION 

Density functional theory @,|] (DFT) is a very pow- 
erful and elegant approach to many-body systems. Its 
foundations are based upon the Hohenberg-Kohn (HK) 
theorem jlj] which states that the total energy of the 
ground-state of a many-particle system in a given exter- 
nal potential can be expressed as a universal functional of 
its electronic charge density. Thus, within a DFT-based 
formalism, one avoids constructing the many-body wave- 
function of the system, and instead computes directly 
ground-state properties like the total energy, magnetiza- 
tion, lattice constant etc. from its charge density. How- 
ever, the exact form of the functional embodied in HK 
is unknown. Therefore approximations are required, the 
most important one being the local-density approxima- 
tion (LDA) developed by Kohn and Sham j| . LDA has 
proven to be extremely successful for the ground-state 
properties of weakly correlated systems. However, for 
strongly correlated systems such as the ones containing f 
electrons, its results are far from satisfactory. There have 
been many improvements proposed to the LDA Q, but 
their applicability depends on the system involved, i.e., 
the present DFT-based approaches are not amenable to 
systematic improvements. 

The alternative approach to the problem of electron 
correlations aims at obtaining the many-particle wave- 
function of the system, by approximately solving the 
corresponding Schrodinger equation within a given one- 
particle/many-particle basis set. Since the Hamiltonian 



of the system under consideration is generally known 
beforehand, there are no uncertainties of principle in- 
volved in this approach. One can systematically improve 
the approach by enlarging the basis set or by includ- 
ing more terms in the expansion of the wavefunction 
of the system. However, such "wavefunction-based" ap- 
proaches scale quite unfavourably with the size of the 
system. Therefore, in their ab initio form, they have 
been applied mainly in quantum chemistry for small to 
medium-sized molecules. In condensed-matter physics, 
wavefunction-based many-body approaches are generally 
used in conjunction with model Hamiltonians, particu- 
larly when the correlations are strong. Usually various 
approximations have to be made [Q and an extension of 
ab initio quantum-chemistry type of methods to infinite 
systems did not seem feasible for a long time. 

Two recent developments may possibly bring ab ini- 
tio wavefunction-based calculations for infinite systems 
within the reach of routine computations: (a) advent of 
more powerful computers (b) use of local excitations in 
electronic structure calculations. The first of the afore- 
said developments is technological in nature, and its im- 
portance is self-explanatory. The second development is 
based on the fact that the correlation hole of an elec- 
tron is a real-space phenomenon. Local or localized or- 
bitals can immediately tell us as to which electrons will 
interact with each other strongly, and which ones weakly. 
Although this has been realized long ago in solid-state 
theory Q as well as quantum chemistry || only recently 
have people started following it in full earnest Re- 
alizing the potential behind the localizcd-orbital-based 
approaches in the electronic structure theory of solids, 
one of us developed the so-called "incremental scheme" 
JlOj l . The incremental scheme is basically an expansion 
of the total corrrelation energy per unit cell of a solid 
written in terms of interactions of increasing complexity 
among the electrons assigned to localized orbitals (Wan- 
nier functions) comprising the solid under consideration. 
Calculations based upon the incremental scheme have 
been performed on a variety of solids JnJ. However, there 
is one possible drawback associated with the computa- 
tions so far, namely, in all correlation calculations the 
infinite crystal was modelled by a finite cluster. For the 
ionic systems the cluster was embedded in a Madelung 
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field, while for the covalent ones the saturation of the 
dangling bonds with hydrogen was the substitute for the 
environment. 

The purpose of the present work is to study electron 
correlation effects in a crystalline insulator within the in- 
cremental scheme, without compromising on the infinite 
nature of the system in any way. Thus by comparing the 
results of the present work with calculations done by the 
previous cluster-type approximation, we can critically ex- 
amine both approaches. Bulk LiH was chosen as a first 
test case for our method. To the best of our knowledge 
this is the first calculation employing a wavefunction- 
based approach that does not truncate an infinite solid 
into a finite cluster. Although the computational efforts 
involved is more than for the cluster-based calculations, it 
is not excessive. Therefore, we believe, that with further 
algorithmic improvements, the present approach can be 
made viable for solids containing many atoms in a unit 
cell. 

The remainder of the paper is organized as follows. In 
section |l| we briefly discuss the theoretical background 
of the present work. In section II] we present the results 
of our calculations performed on bulk LiH. Finally, in 



section IV we present our conclusions. 



II. THEORY 

The starting point of our approach is the representa- 
tion of the restricted Hartree-Fock (RHF) ground state 
of the solid in terms of Wannier functions (WFs). It is 
accomplished by solving the RHF equations of the in- 
finite crystal directly in the Wannier representation, as 
recently described jf2j . The applicability of the approach 
to periodic insulators of reduced dimensions was demon- 
strated in another work of ours |jL3[ , where we used it 
to study the ground states of an LiH chain and trans- 
polyacetylene, at the RHF level. Assuming that the sys- 
tem consists of N unit cells (N —> oo) each of which has 
2n c electrons, its RHF ground state will be described by 
n c doubly occupied WFs per unit cell. If we use Greek 
indices a, (3 etc. to denote the occupied WFs localized 
in a given unit cell, the RHF wavefunction of a crystal 
is a Slater determinant composed of Wannier functions 
{|a(Rj)); a = l,n c ;j = 1,N}. The latter form an or- 
thonormal set 



(1) 



and Wannier functions localized in different unit cells are 
translated copies of each other 



|a(R i + R i )>=r(R i )|a(R i )), 



(2) 



where the operator T(R^) represents a translation by lat- 
tice vector Hi. LiH has 2 occupied WFs per unit cell. 
Although, in our approach, the WFs are assumed to be 



the linear combinations of Gaussian lobe-type functions 
]i"4j , the theory presented here is quite independent of 
the choice of the basis set. 

The next step is to choose a region of space — hence- 
forth to be called the "correlation region" (C) — in which 
the correlation effects will be computed. The infinite 
region of space lying outside the correlation region will 
be called the "environment" {£) . In the correlation cal- 
culations, we consider virtual excitations from the occu- 
pied WFs of the correlation region only, keeping the elec- 
trons of the environment "frozen" to the occupied (RHF) 
space. Because of the frozen character of the environment 
WFs, one can immediately sum up their contribution to 
obtain an effective one-electron potential called the "en- 
vironment potential" 

K7= E (2( P «(R,)i^k«(R,)) 
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where p and q are two arbitrary basis functions and the 
factor of two in the first term is due to the spin sum- 
mation. U^ v is evaluated using the Ewald summation 
technique [^lj, and one is left with an effective Hamilto- 
nian for the electrons located in the region C. The physi- 
cal interpretation of U^ Y is obvious — it represents the 
influence of the electrons of £ on those of C. 

The correlation contribution to the total energy per 
unit cell is computed using the aforesaid incremental ex- 
pansion fnl] 



A ^ = E^ + ^E A ^ + ^ E A 



(4) 



where the summation over i involves Wannier functions 
located in the reference cell, while those over j and k 
include all the Wannier functions of the crystal. The 
6i ("one-body" increment) are computed by considering 
virtual excitations only from the i-th Wannier function, 
freezing the rest of the solid at the HF level. The "two- 
body" increment Aey is defined as Aey = ey — — ej 
where e.y is the correlation energy of the system ob- 
tained by correlating two distinct Wannier functions i 
and j. Thus Ae^ represents the correlation contribution 
of electrons localized on two "bodies" i and j. Simi- 
larly higher-order increments are defined. For ionic sys- 
tems such as LiH, generally one- and two-body incre- 
ments account for the bulk of correlation contribution — 
three-body increments being negligible. The incremen- 
tal expansion of Eq. (Q) is completely independent of the 
method of computation. Since LiH is a small system, 
for the present calculations of increments it was possible 
to use the full configuration-interaction (CI) approach 
which in the physics literature is known as the exact di- 
agonalization method. For larger systems, however, full- 
CI will not be feasible. In that case one can resort to any 
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of the other available size-extensive correlation schemes 

We finally comment upon the construction of the space 
of virtual orbitals used to compute the correlation effects. 
We find it absolutely essential that the virtual orbitals be 
also localized. For this purpose we adopted the approach 
suggested by Pulay Q. It consists of first orthogonaliz- 
ing the basis functions localized in the reference cell — 
using projection operators — to the occupied Wannicr 
functions of the reference cell, and to those located in up 
to third-nearest neighbor unit cells. This yields a set of 
virtual functions orthogonal to the occupied space, but 
not orthogonal amongst each other. Orthogonality in the 
virtual space is finally achieved by a symmetric orthog- 
onalization procedure which does not destroy the local 
character of the functions. In the end, we obtain virtual 
Wannicr functions which are somewhat delocalized due 
to their orthogonalization tails, however, quantitatively 
the derealization is not very significant. Since occupied 
and virtual orbitals in our approach form an orthonormal 
set, almost the whole machinery of electron correlation 
treatments developed during the last decades in ab initio 
quantum chemistry for atoms and molecules can also be 
applied to treat polymers or solids. 

III. CALCULATIONS AND RESULTS 

Among the recent studies of bulk LiH, the most note- 
worthy are HF level studies of Dovesi et al. |l5|] and 
of Bella'iche and Levy |l6|, and DFT-based studies of 
Bellaiche and coworkers ]17| , |l8| |. For the sake of com- 
parison, in our study we used the Gaussian lobe repre- 
sentation of the polarized contracted basis sets used by 
Dovesi et al. [|15|, i.e., [3s, Ip] on hydrogen and [2s, lp] 
on lithium. The experimental fee geometry of the sys- 
tem was assumed with hydrogen on the (0, 0, 0) position 
and Li on (0,0, a/2) position, a being the lattice con- 
stant. First all-electron Wannicr HF calculations were 
performed at fifteen lattice constants in the range 4.0 — 
4.225 A. The results provided us with the optimized lat- 
tice constant and cohesive energy of LiH at the HF level. 
From the orbitals of these calculations the appropriate 
environment potentials were constructed (cf. Eq.([|)), 
which were subsequently used to perform incremental- 
scheme-based correlated calculations, employing the full- 
CI code of Knowles and Handy p9[ . At the correlated 
level we resorted to the frozen-core approximation by in- 
cluding the Li Wannier functions of region C also in J7 env . 
Since Li Is type orbital is highly localized and core-like, 
this approximation is quite reasonable and is analogous 
in spirit to the pseudopotential calculations performed 
routinely in condensed-matter physics. 

In table || we present the results of incremental cal- 
culations performed for the infinite solid for the lattice 
constant of 4.07 A. The increments reported there include 



the one-body increment corresponding to the correlation 
energy of the polarized hydrogen anion of the reference 
unit cell and the two-body increments involving the in- 
teraction of the reference hydrogen anion with those sit- 
uated in up to third-nearest neighbor unit cells. Note the 
rapid decrease of the two-body increments with increas- 
ing distance between the ions which is of van der Waals 
type. The reason behind this can be intuitively under- 
stood as follows. When, for the purpose of the compu- 
tation of correlation effects, atoms (or ions) are virtually 
excited, they develop local (virtual) electric dipole mo- 
ments. Therefore, the two-body increment corresponds 
to the interaction of two localized (virtual) dipole mo- 
ments, the interaction between which is nothing but the 
van der Waals interaction. A noteworthy point is that 
localized orbitals are absolutely essential for this argu- 
ment to hold and, within a Bloch orbital approach, such 
a simplification is not possible. For the sake of compari- 
son, in table | we also present the results of calculations 
performed on a finite cluster modelling the bulk LiH, us- 
ing the same basis set and geometry. The finite cluster 
calculations were performed with MOLPRO program pc| ] 
and include only the most important increments namely 
one-body and the nearest-neighbor two-body increments. 
The maximum disagreement of w 1.6 x 10 -3 Hartrees 
between the two calculations is very encouraging, and 
essentially defines the error bars of the numerics. 

The equilibrium lattice constant, cohesive energy and 
the bulk modulus obtained from our calculations — both 
at the HF and correlated level — are presented in table 
H together with results of other authors pl]-p7fl. Be- 
fore accounting for the zero-point vibrations of the nu- 
clei, the DFT-based calculations systematically underes- 
timate the lattice constant. Since the hydrogen nucleus 
is very light, one would expect the energy of its zero- 
point vibrations to make non-negligible contributions to 
the total energy. Indeed, Roma et al. j2^[ applying den- 
sity functional perturbation theory (DFPT) found that 
the inclusion of the zero-point energy increases the equi- 
librium lattice constant of LiH by 0.08 A. Bellaiche et 
al. [[HI , by adding the lattice-constant-dependent zero- 
point energies obtained by Roma et al. [^7| , to their pre- 
viously derived DFT energies (l7) were able to obtain 
near-perfect agreement with the experimental lattice con- 
stant in case of the Wigner expression for the correlation 
energy, whereas their result is still 0.06 A too small for 
the more modern Perdew-Zunger parametrization of the 
Ceperly- Alder quantum Monte Carlo correlation energy. 
For the wavefunction-based approaches table [n] shows 
that the HF approximation yields a too large lattice con- 
stant due to an overestimation of the ionic radii, par- 
ticularly for the anion. Once electron correlation effects 
are taken into account, we get an agreement of better 
than 0.01 A with the experimental value of the lattice 
constant, without accounting for zero-point effects. 

The experimentally measured cohesive energy of LiH, 
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corresponding to its dissociation into Li + and H~, is 
0.346 a.u. p^| . This, when corrected by the difference 
of the electron affinity of H and the ionization potential 
of Li, leads to a cohesive energy of 0.176 a.u. with respect 
to the separation into free Li and H atoms. In the present 
case, at the correlated level, the Li core is held frozen. 
Thus the theoretical cohesive energies both at the HF 
and the correlated level are obtained by subtracting from 
the total energy per unit cell, the HF limits of the ener- 



gies of free Li (7.4327 a.u. |28|) and H (0.5 a.u.) atoms. 
From table || it is obvious that for the HF cohesive en- 
ergy there is very good agreement among the results of 
all authors. The difference of 0.0006 a.u. between our 
results and the one of Dovesi et al. Jl5|] is possibly due 
to our use of Gaussian lobe functions instead of Carte- 
sian Gaussians. By subtracting our HF cohesive energy 
(0.1294 a.u.) from the experimental one (0.1760 a.u.), 
one gets « 0.0466 a.u. as the "experimental" value of 
the correlation contribution to the cohesive energy per 
unit cell. The correlated cohesive energy that we obtain 
in the present work is 0.1644 a.u. and corresponds to 
93.4 % of the experimental value. The calculated corre- 
sponding correlation contribution of 0.0350 a.u. is 75.1 
% of the experimental value. 

In order to account for a higher percentage of the cor- 
relation contribution we would need to apply larger basis 
sets including additional diffuse and also higher angular 
momentum functions. Unfortunately, since our program 
is still under development, such large-scale calculations 
are not possible at present. Nevertheless, one can esti- 
mate the effect of the basis sets incompleteness for the 
time being from finite cluster calculations using extended 
basis sets |2^]. We obtained a value of -0.035486 a.u. for 
the one-body contribution and -0.001544 a.u. for the 
nearest-neighbor two-body contributions. Since we are 
interested in the values of the corresponding contribu- 
tions obtained with the extended basis sets for the infi- 
nite solid within the Wannier function approach, we scale 
these numbers according to the ratios of the correspond- 
ing numbers in the two columns of table [| and obtain the 
estimates -0.033652 a.u. and -0.001230 a.u., respectively. 
By comparing these values to the corresponding numbers 
obtained with the smaller basis set for the bulk case, we 
obtain corrections of -0.004335 a.u. and -0.000553 a.u. 
for the one-body and two-body term, respectively, lead- 
ing to a total correlation energy gain of -0.0110 a.u. Our 
so-corrected estimate of the correlation contribution to 
the cohesive energy per unit cell of 0.0460 a.u. (0.0350 
a.u. + 0.0110 a.u.) is in almost perfect agreement with 
the corresponding experimental number 0.0466 a.u.). 
We want to point out here, that further technical im- 
provements in our implementation of the method will 
make calculations with extended basis sets for the in- 
finite system feasible and therefore extrapolations from 
finite cluster calculations unnecessary The fact that an 
extension of the basis sets leads to a calculated value 



essentially in agreement with the experimental result un- 
derlines the potential of the proposed method, i.e., the 
possibility of systematic improvement. 

The results of our Wannier-function-based calcula- 
tions, using the smaller basis set of Dovesi et al. JL5| , 
show that in a sufficiently wide range of values around 
the equilibrium lattice constant (15 points in the in- 
terval 4.0 A < a < 4.225 A) to high accuracy (cor- 
relation coefficient 0.999983) the correlation energy de- 
pends linearly on the lattice constant a (E corr (a.u.) — 
-0.034030 + 0.002935 * (a(l) - 4.06)). This is in agree- 
ment with an assumption made in previous work on NiO 
by Doll et al. |ll). Such a "linear law" is important be- 
cause it allows one to perform correlated calculations only 
for two values of the lattice constant near its equilibrium 
value and to deduce the other values by interpolation. 
Assuming also a linear behavior of the correlation contri- 
butions obtained with the extended basis sets mentioned 
above |p9| , we derive an estimated equilibrium lattice 
constant of 4.054 A compared to the actually calculated 
value of 4.067 A and the experimental result of 4.060 A. 

Finally, we discuss the influence of correlations on the 
bulk modulus of LiH. The experimentally measured val- 
ues reported in literature ]23]-|25]| are actually obtained 
by fitting a pressure- volume (p(V)) relationship to ex- 
perimental data points. The values for 7 LiD are only 
about 1 % to 3 % larger than for 7 LiH |2|,|27|] and some- 
times a common value is reported p5| ]. Depending on 
the chosen ansatz for the p(V) equation of state (Birch, 
Murnaghan, Vinet) the values reported for 300 K range 
from 31.1 GPa to 36.2 GPa, with error bars of 1 % to 
6 %. In addition to this uncertainty, the values have to 
be corrected to K, which roughly corresponds to an 
increase by 6 % |24|,|2(|, i.e., a range from 33.0 GPa to 
38.4 GPa. A direct comparison of our calculated values, 
which do not rely on any fitting of an equation of state 
and are for K, and the effective experimental parame- 
ters is to be viewed with care. From table || it is obvious 
that our HF value of 33.4 GPa is in within the interval of 
experimental values, whereas our correlated (estimated) 
value of 36.5 (37.9) GPa appears to be slightly too large. 
Roma et al. |27J observed a reduction of the bulk modu- 
lus by 4 GPa (« 10 %) due to zero-point effects. Using 
their zero-point corrections Bella'iche et al. |l£| were able 
to obtain close agreement with the experimental values 
of the bulk modulus in case of the Wigner correlation 
energy expression (33.9 GPa), whereas for the Perdew- 
Zunger parametrization also a relatively large value was 
obtained (36.1 GPa). We note that by adding the correc- 
tions calculated by Roma et al. |27j] to our energy-volume 
data, we end up with correlated results which are not in 
good agreement with the experimental values, i.e., a lat- 
tice constant of 4.143 A and a bulk modulus of 30.7 GPa. 
Since our uncorrected values appear to be stable with re- 
spect to systematic improvements of the calculation and 
are also quite close to the experimental values, we believe 
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that the zero-point corrections obtained by Roma et al. 
may be slightly too large, despite the fact that previous 
estimates were even significantly larger pCfl . 



IV. CONCLUSIONS AND FUTURE DIRECTIONS 



and by the improvements in the algorithms of electronic 
structure theory. However, in order to be more specific 
about the shapes and sizes of the systems which can be 
studied using the present approach, we believe that more 
experience with it is essential. Work along these direc- 
tions is presently underway in our group. 



In conclusion, a new wavefunction-based approach has 
been presented, which allows the ab initio determination 
of electron correlation effects in a crystalline insulator, 
without abandoning its infinite character. The approach 
merges ideas from quantum chemistry and solid-state 
physics, and has the benefit of being amenable to sys- 
tematic improvements of results. It has been shown to be 
viable by a first successful application to an ionic crys- 
tal, namely bulk LiH. The applicability of the present 
approach to covalent insulators and semiconductors such 
as diamond, Si, GaAs etc. is essentially straightforward. 
The only trivial modification in those cases will be that 
the occupied (and virtual) single-particle states will con- 
sist of bond-centered Wannier functions, instead of the 
atom-centered Wannier functions encountered in case of 
ionic systems. As a matter of fact, the viability of the 
present local approach to electron correlations in cova- 
lent systems has already been demonstrated by Paulus et 
al. [Tl| — albeit within a finite-cluster approximation — in 
their study of Group IV semiconductors. Its generaliza- 
tion within the infinite-crystal method presented in this 
work will be the subject of a future investigation. 

As far as the treatment of excited states of crystals 
is concerned, the quasi-particle excited states such as an 
electron (hole) in the conduction (valence) band can be 
treated straightforwardly using the present approach, as 
these excited states can be regarded as the ground states 
of N + l (N — l) electrons, where N is the total number of 
electrons in the unperturbed system. Of course, the elec- 
tron correlation effects in these cases will be more com- 
plicated as compared to the V-particle ground state, ow- 
ing to the polarization cloud of the extra electron (hole) 
in the system. Applicability of this approach to quasi- 
particle excitations was demonstrated by Grafenstein et 
al. [ jBlp^l , who computed the correlated dispersion rela- 
tions of holes in the valence bands of Group IV semicon- 
ductors. Its generalization to the infinite crystal, as well 
as extension to the conduction-band dispersion relations, 
will also be investigated in future. 

As far as the general computational viability of our 
approach is concerned, undoubtedly, at present, it is 
more demanding than the DFT-based approaches. How- 
ever, we believe that by utilizing modern order-iV type 
methods of electronic-structure theory, it can be made 
competitive with such methods. Conceptually speaking, 
there is nothing in our approach that restricts its appli- 
cability to systems containing unit cells of certain shapes, 
sizes, or number of electrons. The frontiers in these areas 
will be decided by the progress in the computer hardware, 
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TABLE I. Various increments to the correlation en- 
ergy in Hartrees of bulk LiH computed by the Wan- 
nier-function-based approach presented in this work. For 
comparison, the results obtained using the method of Doll 
et al. jll| are presented, for the most important increments, 
under the heading of Finite Cluster. NN stands for nearest 
neighbors. All the calculations were performed at the lattice 
constant of 4.07 A. 



Correlation Increment Wannier Function Finite Cluster 

one-body -0.029317 -0.030915 

two-body (INN) -0.000677 -0.000850 

two-body (2NN) -0.000113 

two-body (3NN) -0.000023 — 



TABLE II. Equilibrium values of various quantities for the 
bulk LiH obtained in this work as compared to those of other 
authors, and the experiment. Quantities under the column 
DFT refer to slightly different exchange-correlation function- 
als used (see Ref. Ilq] for details). The theoretical results 
presented below do not include the effect of zero-point vibra- 
tions, unless specified otherwise. Lattice constants are in A, 
energies in a.u. (i.e., Hartree) and the bulk moduli are in 
GPa. 

Quantity HF° HF b HF C Correlated" DFT d Experiment 
Lattice Constant 4.106 4.102 4.110 4.067 3.91-3.98 4.060 / 

3.99-4.05 e 

Cohesive Energy 0.1294 0.1302 0.1305 0.1644 0.1760 9 
Bulk Modulus 33.4 34.1 35.0 36.5 33. 9-36. l e 33.0-38.4 h 

31.3-36.2* 



This work 
151 



b Ref. 
c Ref. 



d 



Ref. 



only all-electron results are considered. 
e Includes the effect of zero-point motion. 
f Extrapolated T=0 value [glj 

9 The cohesive energy reported here is the "atomic" co- 
hesive energy obtained by correcting the "ionic" cohesive 
energy reported in ref. p2| ]. See text for an explanation. 
h estimated for K by scaling 300 K values with 1.06 

1 300K psW23] . 
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